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ONESAM, A COMPUTER PROGRAM FOR NONPARAMETRIC 
DATA ANALYSIS AND DENSITY QUANTILE ESTIMATION 


by 


Emanuel Parzen and Scott Anderson 
1. Introduction 


ONESAM (one sample) is a FORTRAN mainline program that implements 
the one sample non-parametric data analysis and density quantile estimators 
proposed by Parzen (1979). There are 4 stages to a quantile function and 
density-quantile function analysis of a single data sample. The program: 
Stage 1. Accepts grouped or individual (ungrouped) data, and forms non- 
parametric raw estimates of the quantile function; Stage 2. Forms a raw 
density-quantile function and sample spacings. 

If the user specifies a null hypothesis probability law, the program 
checks the goodness of fit of the data to that probability law, and forms 
non-parametric smooth estimates of the density-quantile function in Stages 
3 and 4. 

Section 2 describes the stages of analysis in more detail. Section 3 
is a listing of the general operating characteristics and requirements of 


the ONESAM program. Section 4 details the procedure to use ONESAM including 


the options available to users. Section 5 lists example IBM JCL setups. 


2. Comments on the Stages of Analysis 


Stage 1: User inputs data, outputs sample quantile function. 


A. If the input is ungrouped data, the user must supply: the tape 
: number (specified by DATA), upon which the data resides (if other than 


unit 5 which is the tape the data is usually on when read in from cards); 


the format of the data; and the data itself. The output of stage 1 con- 


sists of: various descriptive statistics; the order statistics by quar- 


tiles; the empirical quantile function, Q(u), defined by Q(u) = Xien 
a 


for u = cs >» J = 1, ..., nz; and a quantile box plot. 


Options: 


1. Data can be standardized to take values in the unit (0,1) 
interval by specifying IOPT2 = 1. The output (which includes 
descriptive statistics, order statistics by quartiles, and Q) 
is based on the standardized data. When this option is used, 
future stages of analysis will use standardized data for 
analysis. 
Compute Q(u) using linear interpolation, for u = jh, j = 1, -.., NQ 
where NQ = [1/h] by specifying FINC = h. When this option is used, 
future stages will use Q computed at points h, 2h, ..., l-h. 


B. If the user inputs grouped data by specifying IOP1 = 1, the user 
must also supply: the lower limit of the first interval (XMIN); the interval 
width (XINT); and, if desired, the total number of observations (NOBS). The 
user must input the numbers of observations in each interval in Format 
(20F4.0) and supply a value, h (FINC), to be used in computing Q. As with 
the raw data, the user must supply the number of the unit (DATA) that the 


data is to be read from (if other than unit 5). 


Stage 2: Raw density-quantile function £Q(u) and sample spacings. 
In Stage 2, the function f£Q(u), defined by 


2h 


£Q(jh) = ———*¥________ 
Q((4+1)h)-Q((4-1)h 


for u = jh, j = 1, ..., NQ, is computed. Also computed in Stage 2 is 
q(u) = 1/£Q(u) . Plots of these functions are given. 

Note that if Q is not computed via linear interpolation from grouped 
data, but rather is computed from ungrouped data, the definition of 4(u) is 


such that 


a(u) = n{Q((j+.5)/n) - Q((j-.5)/n)} = MX (549) a Xy en! 


6h Sly. 24. sey Red 


Stages 1 and 2 constitute the non-parametric raw analysis of the 
data. If the user supplies a .ull-hypothesis probability law, one can 
proceed to Stages 3 and 4 which constitute a non-parametric smooth 


analysis of the data and also a goodness of fit test for the null hypothesis. 


Stage 3: Raw D(u) and o(v). 


The user must specify a density-quantile function £5% for the null 
hypothesis Hy? Qz=ut Q- This is done using indicator variables. If 
one wishes to test against density "j", then set DQH(j) = 1. The possible 
densities one can use are supplied in Table I. 


The program computes: (1) the raw transformation density, 


d(u) = f9Q)(u)a(u)/dq , 


also called the weighted spacings; and (2) the raw transformation distribution 
function, D(u) called the cumulative weighted spacings. Tests for Hy 

could be obtained by using the facts that under Hy the d's should be dis- 
tributed as exponential with mean 1, and Yn {D(u) - u}is asymptotically a 


a Brownian Bridge stochastic process B(u), or a linear transformation of 


B(u). Py is an estimate of 0, defined by 


a Ii a i i hi a aid ac Ra ihe ll at cs 


| oy 8 if £9 Q% (u) a (u) du : 


The program then computes the Fourier transform 


$(v) = iF d(u) exp(2nivu)du for v = 0, 1, ..., 


ahs 


by specifying ORDER = m (the program chooses order 5 if no other order is 
specified). By comparing the square modulus of the phi's, lacvy) 2 » which 
are plotted, to a suitable threshold value (for example, 2*NOBS), one could 
use these values as evidence to accept or reject Hp: 

The real part of the phi's are submitted as @'(m) , $'(m-1), ..., 
o'(1), «.., o'(m-1), $'(m) to obtain the shifted S array. ‘Two separate 
shifted S array's are computed; first where $'(v) = Re$(v) and second 


where $'(v) = Re(-1)"$(v) for v= 1, 2, ...,m. 


Stage 4: Smooth D(u) and £Q(u). 


Autoregressive smoothed density estimators £Q and smoothed trans- 
formation distribution functions D(u) are computed for successive orders 
1, 2, ..., m. Plots of D vs D are given for all orders to help select ° 
the order that gives the best fit to the data. In addition, 02 ‘ log, (92) ’ 
and the values of CAT and AIC criteria are printed by AR order. If CAT 
(or AIC) selects order 0, one considers this as evidence to not reject 
the null hypothesis. 


Alternative analyses: Stages 3 and 4 can be repeated for several 


possible hypothesized densities simply by setting any other DQH(4) = 1. 


TABLE f 


List of Density-Quantile Functions 


0% 
Normal 

Exponential 
Logistic 

Double Exponential 
Uniform Reciprocal 


Cauchy 


Extreme Value 


.Log Normal 
Pareto 
Weibull 


Half Logistic 


Table II 


A List 


Parameter Value Status Default 


SIZE Real > 0 Optional 0.0 


Integer Optional 


NBATCH 


Integer > 1 Required 


Remarks 
Scale factor for CALCOMP plots. 
If equal to 0.5 the plots are 


half sized. 


Plotting option -- equal to 0 if 
only printer plots are desired. 

Equal to 1 if only CALCOMP plots 
are desired. Not equal to 0 or 1 


if both CALCOMP and printer plots 


are desired. 


Number of data sets to be processed. 
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3. Operating Characteristics 


The ONESAM program is written completely in FORTRAN and requires ap- 
proximately 192 K bytes of storage. The program reads on 1/0 unit 5 and 
writes on 1/0 unit 6. Additionally, I/0 unit 9 is required for scratch 
file purposes and the user may request that the program fetch its data 
set off of another (not 5, 6, 9) I/0 unit. 

Initially, the ONESAM package is set up to require the California 
Computer Products Plotting System (CALCOMP). The program is easily con- 
verted to systems without the CALCOMP plotter by: 

(1) deleting all sections of code enclosed by cards of the form 


C$$$$$$ followed by C 


(2) removing the subroutines GPLOT, LABEL, MXMNSC, and QBXPLT. 


ONESAM requires the inverse normal IMSL subroutine MDNRIS. Systems 
without the IMSL package may substitute an inverse normal subroutine of 
the form 


subroutine MPNRIS (P, Y, IER) 


where P is the input probability value, Y is the output inverse normal 
value, and IER is a dummy error indicator variable. A listing of the IMSL 
subroutine MDNRIS is included in this report on page 75. 
4. Procedures 

This section details the procedure one should use to complete a data 
analysis using the ONESAM package. The most critical feature of the ONESAM 
package is to understand completely the variables listed in Tables II and 
III. These variables contain almost all the information and options used 
in the analysis. 


The procedures are most easily explained by use of the following typical 


JCL setup: 


ee eg Se as Bre ee oe ee ee ee SM Ee | ORR Se ae Pome g ee Pee 


See ae ge oe 


1. 7 PRYRR JOB OOUOG, AAAX, $30,005 NN). “AAAAAN’ -NSELEVEL® (1.0) 

2. /®JOBPARA R219? .K=0 

3. /*LEVEL 0 

4. //JOBLIB DD DSW=VYL.SU.RGD. JOBLIB,DISP=OLD 

3. 77STEP1 EXEC PGR=ONSPLP..REGION=192K 

6. //PLOTLOG BD SYSOUT=A 

7. 7/VECTRI «= - DD DSN=SAVECTRILUNIT=SYSDA.SPACE=(TRK. (5.17) DISP=( PASS) 
8. //VECTR2 BD DSN=SEVECTRZ.UNIT=SYSDA.SPACE=(CYL. (1,17) .DISP=( .PASS) 
9. /7PLOTPARK DD + 


1a. SPLOT NSGLVL=1 .XNIN]-.5.YMIN=-.5.XMAX=810.0 SENT 
\\. //ETOSFOO1 DD + 

12. SA NBATCHEII1 .PPLT=1.3END 

13. BATA SET RANE 

14, RPARMS N=111, SEND 


\5. FORMAT OF DATA EG (14F5.1) 

16. //¥T06FOO1 DD SYSOUT=A 

\7. //FTOFFOO! BU UNIT=SYSBA.SPACE=(TRK.50). 

{8. 77 DCB=(RECFMSFB.LRECL=80 ,BLKSIZE=3120) 

19. JAETOVFOOT DD DSN=WYL XX ARK AXXXX(AAXAX) 
20. /? DISP=SHR.LABEL=(..,1N) 

2\. J/STEP2 EXEC PGH=VTPLOT .REGION=1 92K .COND=(0.8E) 
22, 7/STEPLIB BD DSN=SYS1 VTEC. PLOTLIB.BISP=SHR 
23. 7/PLOTLOG DD SYSOUT=A 

24. //VECTRI «DD DSR=RAVECTRI .DISP= (OLD. DELETE) 
25. FINECTR2 «DD DSN=RSVECTR2. DISP=(OLD DELETE) 
26. /ISNSVECTR DD SYSOUT=P .DCB=LRECL=1 26 

27. /*END 


Cards 6 - 10 and 21 - 26 are only required if the CALCOMP plotting is 
requested, i.e., PPLT #0. Cards 12 and 14 utilize the namelist feature 
of the FORTRAN language and contain the user specified options of Tables II 
and III respectively. Both of these cards are required even if no options 
are requested. Card 12 must begin with b&Ab and end with ,&END, 

(Note: b denotes a blank space). Card 14 must begin with b&PARMSb 

and end with ,&END. User options are specified as illustrated. Card 15 
is required only if the user inputs ungrouped data. Data format is assumed 
(20F4.0) if grouped data. Card 11 specifies the input unit. Card 16 
specifies the output unit. Cards 17 and 18 specify the scratch unit. Cards 
19 and 20 specify the optional unit on which the user has his data file 
located. If the user wishes to include his data set within his JCL deck 


he would list his data between cards 15 and 16. 


Cards 12 - 15 are repeated NBATCH times. 


=ji= 


5. Sample JCL 
Sample JCL are given for the following types of jobs: 


i A. Ungrouped data on external file (WYLBUR) using I/0 unit 1; 
B. Ungrouped data on cards, with three sets of data; 


C. Grouped data on cards and requesting CALCOMP plotting. 


-12- 


A/SCOTT JOB (RS79.0028 .530.005.5U).ONSPLP’ .NSGLEVEL=(1.0) 
/*JOBPARK R=192.Ka0 

7*LEVEL 0 

/7JOVLIB BD BSNVYL.SU ROD. JOBLIB .DISP=0L8 
/7STEP1 «= EXEC PGR=ORSPLP REG TOR|1 92K 
//FTOSFOO1 BD + 

SA SEND 
BUFFALO SROVFALL BATA (1910-1972) 

SPARNS N=63,DGH(1)=1 DATA=1 .FINCSO.05. SEND 
(5F 10.0) 

/7FTOSFOO1 DB SYSOUT=A 

7/FTO9FOO! DB URTT=SYSDA .SPACE=( TRK.50), 

‘7 BCB=(RECFN=FB.LRECL*80. BLKSIZE=31 20) 
//FTOVFOO1 DD DSR=UYL AU. ENP QUANT ( BFSHO). 
// DISP=SHR.LABEL=(,..18) 

/END 


-13- 


piers i 
a eee iat Ba 


1.  7/SCOTT JOB (R579 ,002B,530.005,SW),“ONSPLP’ ,NSGLEVEL=(1,0) 
2: /SJOBPARN RE192.KEO : 
3. /#LEVEL 0 
a. 7/SOBLIB DD DSNeVYL.SU.ROD. JOBLIB, DISP=OLD 
5. VISTEP1 EXEC PGM=ONSPLP.REGION=192K 
7 &. //FTOSFOO1 DD * 
2. SA NBATCH=3, BEND 
8. _- TIPPETT’S WARP BREAKS AL 
9. SPARNS N=9,DOR(1)=1 .DOH(2)=1 , ZEND : 
10, (10FS.0) 
: iW. 26. 30. 54. 25. 70. 52. Sl, 26. 67. 
‘ 12. _ TIPPETT’S WARP BREAKS AN 
13. SPARMS N=?.DQH(1)=1,DON(2)=1 SEND 
14. (10FS.0) 
18. 18. 21. 29. 17. 12. 18. 35. 30. 36. 


1d. TIPPETT’S WARP BREAKS AH 
\7. SPARNS N=9.DQH(1)=1 .DGH(2)=1 SEND 

: 8. (10F5.0) 
V9. 36. 21. 24. 18. 10. 43. 28. 15. 26. 
20. 7/FT04FOO1 DD SYSOUT#A 
2\. //FTQ9FO01 DD UNIT2SYSBA,SPACE=(TRK,50), 
22. 77 DCB2(RECFM=FB.LRECL#80.8LKSIZE23120) 
23. /8ENB 


| JISCOTT JOB (R579 .0028.530,005.SU).’ONSPLP’ .ASGLEVEL= (1.0) 


/*JOBPARH R=196.K20 
ALEVEL 0 
//JOBLIB DS DSN=VYL.SU.ROD. JOBLIB.DISP=0LD 
/7STEP1 EXEC PGM=ONSPLP.REGION|196K 
7/PLOTLOG DD SYSOUT=A 
7INECTRI DB DSMERAVECTRI UNL T=SYSBA.SPACE=(TRK,(S.1)) .BISP=( PASS) 
ZIVECTR2 BD DSW=SRVECTR2.UNIT=SYSDA SPACE=(CYL. (1.12) DISP=( PASS) 
/7PLOTPARN BD # 

SPLOT MSGLVL=1 .XNINS-.5,YHIN®-.5.XNAX=810.0 SEND 

//ETOSFOO1 DD * 

3A PPLT=1.2END 
BLOOD PRESSURE DATA--NOW-USERS--AGE 17-24 

BPARMS N=16,DQH(1}=1. JOPTI=1 XMIN=OS. .XINT=S. .PINC!.05. SEND 

0 1 3 1 t1 18 2013 10 8 4 3 2 0°90 «0 | 

/IETOSEOO1 DD SYSOUT=A 
7/FTOVFQO1 DD UNITSSYSDA SPACE=(TRK.50). 
/? DCB=(RECF NFB, LRECL #80. BLKSIZE=31 20) . 
//STEP2 EXEC PGN=VTPLOT .REGION=1 96K. COND=(0.NE) | 
//STEPLIB DD DSN=SYS1 .VIEC.PLOTLIB. DISP=SHR 
H7PLOTLOG DD SYSOUT=A 

/NECTRI DD BSN=SSVECTRI.DISP= (OLD, DELETE) 
Z/VECTR2 BD DSN=S&VECTR2. DISP= (OLD, DELETE) 
VISYSVECTR DD SYSGUT=P.DCB=LRECL=1 24 

/AEND 


-15- 
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